library(sfsmisc)
library(rgl)

rsm <- 1
rsmdata <- read.csv("/home/rumen/cryptsim/analysis/run.simulation.coefficients.models.csv",header=FALSE)

rsmmodel <- rsmdata[rsmdata[,1]==rsm,]
repr <- as.numeric(strsplit(as.character(rsmmodel[1,3])," ")[[1]])
surv <- strsplit(as.character(rsmmodel[2,3])," ")[[1]]
muta <- strsplit(as.character(rsmmodel[3,3])," ")[[1]]

repr.reduced <- repr[repr!=0]
relative.fitness.coefficients <- NULL 
for(i in 1:length(repr.reduced)){
 combos <- combn(repr.reduced, i)
 combos <- combos + 1
 relative.fitness.coefficients <- c(relative.fitness.coefficients, apply(combos, 2, prod))
} 
reproduction.fitness.classes <- c(1.00,sort(unique(relative.fitness.coefficients)))

# Prepare fitness landscape 4 by 8 = 32 cells
yy <- expand.grid(locus1 = c(0,1),  locus2 = c(0,1))
xx <- expand.grid(locus3 = c(0,1),  locus4 = c(0,1), locus5 = c(0,1))

combos <- digitsBase(0:31,base=2,n=5)
all.loci.combos <- (repr.reduced+1) * combos
all.loci.combos[which(all.loci.combos==0)] <- 1.0
fitness.combos <- apply(all.loci.combos,2,prod)

orrd <- order(fitness.combos)
sorted.combos <- combos[,orrd]
sorted.fitness.combos <- fitness.combos[orrd]
#all.loci.combos <- all.loci.combos[,orrd]

z <- matrix(sorted.fitness.combos,nrow=4,ncol=8,byrow=TRUE) 
rg.ramp <- colorRampPalette(c("darkgreen", "yellow", "red","blue"))
ncolors <- 32
col.arr <- rg.ramp(ncolors)
fitness.col <- cbind(sorted.fitness.combos, col.arr)

f1 <- c(0.3352588 ,-0.940788, -0.05019309   , 0)
f2 <- c( 0.4839025 , 0.126242 , 0.86596853  ,  0)
f3 <- c(-0.8083563, -0.314612 , 0.49757338  ,  0)
f4 <- c(0.0000000,  0.000000,  0.00000000  ,  1)
fov <- rbind(f1,f2,f3,f4)
open3d(userMatrix=fov,windowRect=c(0,0,350,250), viewport=c(0,0,350,250),zoom=1.1)

# Scale down by a factor so that the maximum height is 4, 0.02 is background division rate
wt.div.rate <- 0.02
newz <- (z*wt.div.rate) * (4.0/max(z*wt.div.rate))
zoffset <- 0
shift <- 0.1
for(i in 1:nrow(z)){
for(j in 1:ncol(z)){
      # Make the 4 top points
      ttl <- c(i-1+shift,j-1+shift,newz[i,j])
      ttr <- c(i-1+shift,j-shift,newz[i,j])
      tbl <- c(i-shift,j-1+shift,newz[i,j])
      tbr <- c(i-shift,j-shift,newz[i,j])
      # Make the 4 bottom points 
      btl <- c(i-1+shift,j-1+shift,zoffset)
      btr <- c(i-1+shift,j-shift,zoffset)
      bbl <- c(i-shift,j-1+shift,zoffset)
      bbr <- c(i-shift,j-shift,zoffset)
      # Make 1 quad that is the top
      cover <- rbind(ttl,ttr,tbr,tbl)
      # Make the walls
      w1 <- rbind(btl,btr,ttr,ttl)
      w2 <- rbind(btr,bbr,tbr,ttr)
      w3 <- rbind(bbr,bbl,tbl,tbr)
      w4 <- rbind(bbl,btl,ttl,tbl)
      # Make walls transparent
      quads3d(rbind(w1,w2,w3,w4),col="lightgrey", line.antialias=TRUE)
      quads3d(cover, col=fitness.col[z[i,j]==fitness.col[,1],2],lit=FALSE, line.antialias=TRUE,point.antialias=TRUE)
}
}
a2 <- format(2/4*max(z*wt.div.rate),digits=2)
a4 <- format(max(z*wt.div.rate),digits=2)
axis3d(c("z+-"),tick=FALSE, at=c(0,2,4),labels=c(0,a2,a4),line=1)
lines3d(rbind(c(0,8,2),c(4,8,2),c(4,0,2)),col="black",lit=FALSE)
# make bounding box comprised of 3 walls
# bottom plate
quads3d(rbind(c(0,0,0),c(4,0,0),c(4,8,0),c(0,8,0)),col="black", line.antialias=TRUE,front="lines")
# left plate
quads3d(rbind(c(0,8,0),c(4,8,0),c(4,8,4),c(0,8,4)),col="black", line.antialias=TRUE,front="lines")
# farther plate
quads3d(rbind(c(4,0,0),c(4,0,4),c(4,8,4),c(4,8,0)),col="black", line.antialias=TRUE,front="lines")
mtext3d(text="Genotypes","y--",line=1)
text3d(c(8,1,3),text="Bifurcation rate (b)")
rgl.postscript("fitness_landscape.pdf","pdf")
